\documentclass[9pt,twoside]{osajnl}
%% Please use 11pt if submitting to AOP
% \documentclass[11pt,twocolumn,twoside]{osajnl}

\journal{optica} % Choose journal (ao, aop, josaa, josab, ol, optica, pr)
\usepackage{subfigure}
\usepackage{graphicx}

\usepackage{amsmath}
\usepackage{cleveref}
\DeclareMathOperator{\sinc}{sinc}
\DeclareMathOperator{\somb}{somb}
% See template introduction for guidance on setting shortarticle option
\setboolean{shortarticle}{false}
% true = letter / tutorial
% false = research / review article
% (depending on journal).
\usepackage[justification=centering,
            format=plain]{caption}

\title{Caculation of SFR}

\author[1,2]{Haotian Song}
\author[1,3,4,$\dagger$]{Xiaoyu Nie}
\author[1]{Zhaoyu Zuo}
\affil[1]{Physics Department, Xi'an Jiaotong University, Xi'an, Shaanxi 710049, China}
\affil[2]{School of Physics \& Astronomy, University of Manchester, Manchester M13 9PL, United Kingdom}
\affil[$\dagger$]{nxy1125@tamu.edu}
\setboolean{displaycopyright}{true}
\begin{abstract}
This PDF show the process of caculating SFR and $\delta N$
 
\end{abstract}
\begin{document}
\maketitle
\section{Methods}

In order to analyze ULX populations, we utilized the Binary Population Synthesis (BPS) code \cite{hurley2002evolution} with further updates, which are briefly described as follows. 
% We set-up a grid of initial binary parameters as following.
% \begin{equation}\begin{aligned}
%     M_{1} &: \quad 2 \rightarrow 150 \mathrm{M}_{\odot} \\
%     q &: \quad 0.08 \rightarrow 1 \mathrm{M}_{\odot} \\
%     a &: \quad 3.0 \rightarrow 10^{4} \mathrm{R} \odot
%     \end{aligned}
% \end{equation}
The Initial Mass Function (IMF) developed by \cite{kroupa2001variation} was applied. Although the shapes of IMF differ significantly, it is less subject to the uncertainty due to the similarity of the IMFs for larger masses. The mass of star on zero-age sequence (ZAMS) satisfies the following distribution with initial mass ranging from 2 $M_\odot$ to 150 $M_\odot$.
\begin{equation}
    \xi(m) \propto m^{-\alpha}
\end{equation}
where
\begin{equation}
    \alpha(m)=\left\{\begin{array}{ll}
        +0.3 \pm 0.7 & 0.01 \leq m < 0.08 \\
        +1.3 \pm 0.5 & 0.08 \leq m < 0.50 \\
        +2.3 \pm 0.3 & 0.50 \leq m < 1.00 \\
        +2.3 \pm 0.7 & 1.00 \leq m
           \end{array}\right.
\end{equation}
The distribution of mass ratios q ($q=M_1/M_2$) was uniform between 0.08 and 1 by steps of 0.10 \cite{wiktorowicz2019observed}. The orbital separation is ranged from 3.0 to $10^4 R_\odot$ which facilitates the comparason with other models \cite{Yungelson, Livio & Tutukov
1997}. Here we assume that $ln(a)$ is uniformly distributed. The maxium evolution time is 200Myr in isolation because young galaxies were taken for comparison later. As for supernova kicks, we draw them from a Maxwellian distribution with $\sigma=265km s^-1$ \cite{hobbs2005statistical}. And solar metallicity $Z=0.5Z_\odot$ was adopt in our simulations. 

For sources undergoing Roche Lobe Overflow (RL) mass-transfer (MT), it is estimated that the luminosity of X-ray is \cite{shakura1973black} \cite{poutanen2007supercritically}
\begin{equation}L_{X}=\left\{\begin{array}{ll}
    L_{\mathrm{Edd}}\left(1+\ln \dot{m}_{\mathrm{tr}}\right) & \dot{m}_{\mathrm{tr}}>1 \\
    L_{\mathrm{Edd}} \dot{m}_{\mathrm{tr}} & \dot{m}_{\mathrm{tr}} \leqslant 1
    \end{array}\right.
\end{equation}
where $\dot{m}_{\mathrm{tr}}=\dot{M}_{\mathrm{tr}} / \dot{M}_{\mathrm{Edd}}$ and $\dot{M}_{\mathrm{tr}}$ is in units of $M_\odot/yr$. According to \cite{king2008accretion}, the apparent luminosity is 
\begin{equation}
    L_{app}=L_X/b
\end{equation}
where b is beaming factor defined as $b \stackrel{\text { def }}{=} \Omega / 4 \pi$, and $\Omega$ is the combined solid angle of both beams. \cite{king2009masses} implies
\begin{equation}
    b=\left\{\begin{array}{ll}
    1 & \dot{m}_{\mathrm{tr}} \leqslant 8.5 \\
    \frac{73}{\dot{m}_{\mathrm{tr}}^{2}} & 8.5<\dot{m}_{\mathrm{tr}}
    \end{array}\right.
\end{equation}
whereas for $\dot{m}<8.5$, it is assumed unbeamed. And this formula have been verified in plenty of hyperluminous X-ray sources \cite{king2014hlx,king2016ulxs,king2017pulsing}. In addition, there was an implication \cite{wiktorowicz2017origin} that extremely beamed sources are scarced and transient that lead to difficulty in detection. That means there is no lower limit for the beaming factor.

Furthermore, the wind Roche lobe overflow (WRL)could be taken into acount. The traditional Bondi-Hoyle-Lyttleton (BHL) MT formula \cite{bondi1944mechanism,edgar2004review} is 

\begin{equation}
    \dot{M}_{\mathrm{BHL}}=\pi R_{\mathrm{BHL}}^{2} v_{\mathrm{rel}} \rho
\end{equation}

where $\dot{M}_{\mathrm{BHL}}$ is the mass accretion rate, $v_{\mathrm{rel}}=\sqrt{v_{\beta}^{2}+(v_{\text {orb}}[q /(1+q)])^{2}}$ is relative speed between the wind and compact object. $v_{\text {orb}}$ is orbital speed given by $v_{orb}=2\pi a / P_{orb}$ with $P_{orb}$ being the orbit period. $R_{BHL}=2 G M_{\bullet} / v_{\mathrm{rel}}^{2}$ is the modified accretion radius and $\rho$ is the density of the wind at orbital separation. After assumption of isotropic dilution of the wind, the fraction of wind captured $\mu_{\mathrm{BHL}}=\dot{M}_{\mathrm{BHL}} / \dot{M}_{\star}$ can be:
\begin{equation}
    \mu_{\mathrm{BHL}}=\frac{(1+q) / q^{3}}{\eta(1-f \mathcal{E})^{\beta}\left[1+\left(\eta(1+q)(1-f \mathcal{E})^{\beta} / q\right)\right]^{3 / 2}}
\end{equation}
where $\mathcal{E}$ is the ratio of stellar Roche lobe radius by the orbital separation and only related to q \cite{eggleton1983aproximations}, and $\eta$ is the speed ratio. Another research \cite{el2019wind} proposed that ULXs could remain stable in WRL with a highly beamed wind. Based on realistic acceleration profiles, they estimated the MT by computing the bulk motion of the wind. By assuming non-isotropic dilution of stellar wind, the fraction of wind captured can be significantly higher especially for higher wind speed ratio. And several binary systems have been proved possible for WRL, such as M101 \cite{liu2013puzzling},and P13 \cite{furst2018tale}. We adopted his simulation and assume $\mu$ as a function of $\eta$, stellar filling factor, $\beta$ and mass ratio $q$. The WRL efficiency systematically improved. 

As for the observation of ULXs, the statistics of ULXs in seven galaxies have been published \cite{wolter2018x}. Relevant properties are listed in Table???. For numerical comparison with them, the evolved population of binaries should be in conjunction with a realistic birth rates and metallicity. Here the distribution of each binary system is
 \cite{hurley2002evolution}
 \begin{equation}
    \label[]{eq:hurley}
    \delta r_{j}=\overline{S_b} \Phi\left(\ln M_{1 j}\right) \varphi\left(\ln M_{2 j}\right) \Psi\left(\ln a_{j}\right) \delta \ln M_{1} \delta \ln M_{2} \delta \ln a
\end{equation}
where r is the rate of particular population, and $\overline{S_b}$ is the star forming rate (SFR) in units of number/yr which can be obtained from following formula with SFR in units of $M_\odot/yr$ for binaries massive than 5 $M_\odot$ \cite{grimm2003high}:
\begin{equation}
    \overline{S_b}=(\frac{M_{low}}{5})^{1-\alpha}\frac{(\alpha-2)SFR}{5(\overline{3}+\overline{q})(\alpha-1)}
\end{equation}
where $M_{low}$ is the lower limit initial mass for BPS. For summation of seven galaxies $SFR\approx 42 M_\odot/yr$, $\overline{S_b}=1.5039 yr^{-1}$ The total number N should be the summation of every binary system distribution.
\begin{eqnarray}
    N=\sum \delta N = \sum \delta r \delta T 
\end{eqnarray}
where $\delta T$ is the duration of particular binary evolution stage. 














\begin{equation}
    L_{\mathrm{X}}=\left\{\begin{array}{ll}
    L_{\mathrm{Edd}}\left(1+\ln \dot{m}_{\mathrm{tr}}\right) & \dot{m}_{\mathrm{tr}}>1 \\
    L_{\mathrm{Edd}} \dot{m}_{\mathrm{tr}} & \dot{m}_{\mathrm{tr}} \leqslant 1
    \end{array}\right.
\end{equation}

where $\dot{m}_{tr}=\dot{M}_{tr}/\dot{M}_{Edd}$
\begin{equation}
    L_{edd}=\eta \dot{M}_{Edd} c^2
\end{equation}
whether $\eta$ is 1 or 0.1,

















I am not sure whether SFR includes all stars or binaries only, so I caculate it individually.

$S_b$and $S_s$ are in units of Num/yr.

SFR is in units of $M_\odot/yr$

The main idea of caculation is as following.
\begin{equation}
    \label{eq:1}
    S_b \dot (average of M_{low} and M_{up}) + S_b \dot (average M_2) + S_s \dot(average mass) = SFR
\end{equation}

Equation\ref{eq:1}can be written as following, too.

\begin{equation}
    S_b \dot \int_{M_{low}}^{M_{up}} m_1 \epsilon({m_1})dm_1 + S_b \int_{M_{low}}^{M_{up}} \epsilon({m_1}) (\int_{M_{low}}^{M_1}\frac{m_2}{m_1}  dm_2) dm_1 + S_s \int_{M_{low}}^{M_{up}} m \epsilon({m})dm = SFR
\end{equation}


\begin{equation}
    S_b \dot \int_{M_{low}}^{M_{up}} m_1 \epsilon({m_1})dm_1 + S_b \int_{M_{low}}^{M_{up}} \epsilon({m_1}) m_1 \overline{q}  dm_1 + S_s \int_{M_{low}}^{M_{up}} m \epsilon({m})dm = SFR
\end{equation}

chosen q in my caculation is 0.08, 0.18, 0.28... 0.98, so $\overline{q}=0.535$

If I chose $M_{low}=5M_\odot$, as\cite{grimm2003high} .
\begin{equation}
    SFR=(S_b(1+\overline{q}) + S_s)\int_{M_{low}}^{\infty} m \epsilon(m)dm
\end{equation}

However, the Initial Mass Function differs.
When $m>1M_\odot$, $\epsilon(m)=am^{-\alpha}$

For \cite{kroupa2001variation}, $\alpha=+2.3\pm0.7$

For \cite{kroupa1993distribution} $\alpha=+2.7$

\begin{equation}
    SFR=(S_b (1+\overline{q}) + S_s)\int_{M_{low}}^{\infty} am^{1-\alpha}dm
\end{equation}
\begin{equation}
    SFR=(S_b (1+\overline{q}) + S_s) a\frac{m^{2-\alpha}}{2-\alpha}\bigg|_{M_{low}}^{\infty}
\end{equation}
\begin{equation}
    SFR=(S_b (1+\overline{q}) + S_s) a\frac{M_{low}^{2-\alpha}}{\alpha-2}
\end{equation}

For my BPS caculation, $M_{low}=2M_\odot$, so
\begin{equation}
\overline{S_b}=S_b\frac{\int_2^\infty\epsilon(m)dm}{\int_5^\infty\epsilon(m)dm}=S_b(\frac{2}{5})^{1-\alpha}
\end{equation}

And the relation between $S_b$ and $S_s$ is 
$\frac{2S_b}{S_s+2S_b}=f=0.5$
where f is the fraction of binaries and all stars. If $f=0.5$, then $S_s=2S_b$

And we will get a in the following.

As is shown before,  $\epsilon(m)=am^{-\alpha}$. And  $\epsilon(m)$ should be normalized from $M_{low}$ to $M_{max}$, so a is changed and not equal to the Initial a in the essay.
$$1=\int_{M_{low}}^{M_{max}}am^{-\alpha}dm$$

$$1=a\frac{m^{1-\alpha}}{1-\alpha}\bigg|_{M_{low}}^{\infty}=a\frac{M_{low}^{1-\alpha}}{1-\alpha}$$

$$a=\frac{\alpha-1}{M_{low}^{1-\alpha}}$$

In conclution, $S_b$ can be written as following.
\begin{equation}
    \overline{S_b}=(\frac{2}{5})^{1-\alpha}\frac{(\alpha-2)SFR}{(\overline{3}+\overline{q})M_{low}(\alpha-1)}
\end{equation}
where $\overline{3}$ means whether SFR contains single star. If not, $\overline{3}=1$.

\begin{table}[h!]
\centering
\begin{tabular}{lll}
         & $\alpha=2.3$ & $\alpha=2.7$ \\
Binary   & 4.1560       & 10.6984      \\
All star & 1.8046       & 4.6455
      
\end{tabular}
\end{table}
%\noindent\textbf{Disclosures.} The authors declare no conflicts of interest.
\section{Caculate $\delta N$ }
According to \cite{hurley2002evolution}, the $\delta r$ can be written as:
\begin{equation}
\delta r_{j}=S \Phi\left(\ln M_{1 j}\right) \varphi\left(\ln M_{2 j}\right) \Psi\left(\ln a_{j}\right) \delta \ln M_{1} \delta \ln M_{2} \delta \ln a
\end{equation}
and $\delta N = \delta r \delta T$, so
\begin{equation}
    \delta N = \overline{S_b} M_1 \epsilon({m_1}) \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a \delta T
\end{equation}
And the q and ln(a) are mean distribution, so $k_q\delta q = 1/(q_{number})$,$k_a\delta \ln a=1/a_{number}$
where $q_{number}$ and $a_{number}$ is the cycle index in my caculation. $1/a_{number}$ is the length of step after normalization.
\begin{equation}
    \delta N = \overline{S_b} a M_1^{1-\alpha} \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a \delta T
\end{equation}

\begin{equation}
    \delta N = \overline{S_b} (\alpha-1) (\frac{M_1}{M_low})^{1-\alpha} \frac{\ln(M_{max})-\ln(M_{min})}{M_{number}} \frac{1}{q_{number}} \frac{1}{a_{number}} \delta T
\end{equation}










\section{Adding A to the formula}

The sum of all Success cases should be 1, 
\begin{equation}
    \Sigma Success  = 1,\iiint M_1 \epsilon({m_1}) \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a=1
\end{equation}
but it is not because of failure cases. So A is added.
\begin{equation}
    \Sigma Success + \Sigma Failure  = A \iiint M_1 \epsilon({m_1}) \delta \ln M_1 k_q \delta q \dot k_a \delta \ln a 
\end{equation}
\begin{equation}
    \Sigma Success + \Sigma Failure  = 1 + \Sigma Failure
\end{equation}
\begin{eqnarray}
    A=\frac{1}{ \iiint - \Sigma Failure}=A
\end{eqnarray}

\begin{equation}
    \delta N =(1+\Sigma Failure) \overline{S_b} (\alpha-1) (\frac{M_1}{2})^{1-\alpha} \frac{\ln(M_{max})-\ln(M_{min})}{M_{number}} \frac{1}{q_{number}} \frac{1}{a_{number}} \delta T
\end{equation}

And I got $A=1.0295$
\bibliography{sample}


\end{document}